function [irf,irfd]=girf(sysm,e,horizon,replic1,replic2,cumidx)
% [irf,irfd] = GIRF(SYSM,E,HORIZON,REPLIC1,REPLIC2)
%  calculates generalized impulse response 
%  based on Koop, Pesaran, and Potter (1996, Journal of Econometrics)
%
% INPUT:
%   E       : normalized initial shock
%   HORIZON : length of horizon for responses
%	REPLIC1 : number of current shocks (impluses) to be averaged over
%	REPLIC2 : number of paths of future shocks
%
% OUTPUT
%	IRF		: classical representation of generalized impulse-response
%	IRFD	: generalized impluse repsponse
%
%	IRFD(vidx,hidx,:) gives response distribution (realizations)
%					  for vidx variable at hidx horizon
%   IRF = std(IRFD,0,3);

% Sungbae An
% Updated: 09/06/2005


% check input
if nargin==5
	cumidx=0;
end

% dimension
nz = size(sysm.R,1);
[ny,ns] = size(sysm.ZZ);
nx = ns - nz;

% storage
irfd = zeros(ny,horizon+1,replic1);

% initial shocks
u = randn(length(e),replic1);
u = diag(e)*u;

% initial shock loop
for r = 1:replic1
	
	% storage
	y = zeros(ny,horizon+1);

	% current shock
	w = u(:,r);

	% generate future shocks
	v = randn(nz,replic2,horizon+1);

	% x's are all zero at steady state
	x1 = zeros(nx,1);
	x2 = repmat(x1,1,replic2);

	% initial response
	z1 = sysm.R*w;
	z2 = sysm.R*v(:,:,1);
	x1 = compst(z1,x1,sysm);
	x2 = compst(z2,x2,sysm);

	z1 = repmat(z1,1,replic2);
	x1 = repmat(x1,1,replic2);

	y(:,1) = sysm.ZZ*mean([x2-x1; z2-z1],2);

	% horizon loop
	for n = 2:size(v,3)

		ans1 = sysm.R*v(:,:,n);

		z1 = sysm.T*z1 + ans1;
		z2 = sysm.T*z2 + ans1;

		x1 = compst(z1,x1,sysm);
		x2 = compst(z2,x2,sysm);

		y(:,n) = sysm.ZZ*mean([x2-x1; z2-z1],2);

	end		% n loop

	if sum(cumidx)>0
		% cumulative response
		y(cumidx,:) = cumsum(y(cumidx,:),2);
	end

	irfd(:,:,r) = y;

end % r loop

irf = std(irfd,0,3);

